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A high performance model Stirling engine, in which the heater, regenerator and cooler as a whole is 
formed by hundreds of porous metal sheets, is identified for theoretical analysis to facilitate the future 
scale-up design. The reciprocating flow and heat transfer both in the heat exchanger and in the full 
engine is simulated by a dynamic mesh Computational Fluid Dynamics (CFD) method, and is validated by 
analytical solutions and experimental data. An optimization method is also developed to incorporate the 
entropy generation caused by flow friction and irreversible heat transfer. The results show that relatively 
high indicated power of 33.4 W is obtained, corresponding to a specific power of 1.88 W/cm 3 and a 
thermal efficiency of 43.9%, which are attributable to the extremely small flow friction loss and excellent 
heat transfer characteristics in the regular shaped microchannels, as well as to the compact heat 
exchanger design that significantly reduces the dead volume. Given the same operating conditions, the 
optimized porous-sheets regenerator has a significantly lower total loss of available work while main¬ 
taining even higher thermal effectiveness in comparison with the optimized conventional wire mesh 
regenerator. 

© 2013 Elsevier Ltd. All rights reserved. 


1. Introduction 

Stirling engine has a wide application prospect in the power 
generation field due to its many advantages including adaptability 
to versatile heat sources, high thermal efficiency and environ¬ 
mental friendliness [1 . However, the relatively low specific power 
compared to that of the internal combustion engines is still one of 
the major obstacles hindering its development. So far, extensive 
research has been done on the regenerator, the central and crucial 
component of the Stirling engine in order to improve the engine 
performance [2], 

The traditional wire mesh type regenerator is most popularly 
adopted in Stirling engines due to its huge heat transfer area, high 
convective heat transfer coefficient brought by the cross flow 
around numerous cylindrical shaped wires, and low axial thermal 
conductance. However, there are some inherent disadvantages 
associated with the wire mesh type regenerator [3], such as: (1) the 
numerous cylinders in cross flow produce flow separation, wakes, 
eddies and stagnation zones, resulting in high flow friction and 
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considerable thermal dispersion, a loss mechanism that increases 
apparent axial conduction, damaging power output and engine 
efficiency; (2) the wire screens have some randomness in stacking, 
causing locally non-uniform porosity and flow distribution, which 
might increase axial conduction and damage its thermodynamic 
performance; (3) the mesh wires are subject to the impact of high¬ 
speed high-frequency oscillating flow during operation, so there 
exists the possibility of working loose or fiber breakage, thus 
damaging vital engine components; (4) the wire mesh type 
regenerator also requires long assembly time which tends to in¬ 
crease their cost. 

Theoretically a regenerator with heat transfer surfaces parallel 
to the oscillating flow has a better performance than wire mesh 
type regenerators [4], With the emerging micro-fabrication tech¬ 
niques, properly designed regular-shaped microchannel type 
regenerator can be fabricated to obtain extremely low flow friction 
while maintaining high heat transfer. The main features of the 
regular microchannel type regenerator include: (1) the heat 
transfer surface is smooth; (2) the flow acceleration rates are 
controlled; (3) the flow separation is minimized; (4) the axial 
thermal conduction is reduced by interrupting the axial continuity 
of solid structure, for example, using porous sheets with interme¬ 
diate gaps or clearance. Other advantages include improved 
structural durability, no gas leakage or short-circuit loss owing to 













32 


Z. Li et al. / Energy 64 (2014) 31-43 


tight tolerance, low cost for mass production, thus realizing 
significantly higher comprehensive performance [5], 

Researchers in the US has conducted a series research on the 
regular microchannel type regenerator under the NASA support [6], 
The oscillating-flow rig test showed the highest figures of merit 
ever recorded for any regenerator tested in that rig over its ~ 20 
years of use, demonstrating a shift strongly in the direction of the 
theoretical performance of ideal parallel-plate regenerators. Nu¬ 
merical projections of engine performance using the “SAGE” code 
indicated a performance improvement by 6%—9%. Takizawa et al. 
[7] developed a porous sheet type regenerator with electrically 
etched holes. Performance test in a 3-kW Stirling engine shows that 
the engine performance was improved by about 5%— 10% compared 
to conventional stacked wire mesh regenerator. A series of engine 
tests were also done by Matsuguchi et al. [8] to optimize the 
geometrical parameters of the porous-sheets regenerator. 

Nam et al. [9] developed a parallel wire type regenerator. The 
axial conduction loss is alleviated by wire segmentation, but the 
number of segmentation is limited for the parallel wires. 

Commonly used traditional methods for analyzing and designing 
Stirling engines include the first, the second and the third order 
analysis [10]. Recently, Cheng and Yang [11] developed a lumped- 
mass analytical model to determine the temperature variations in 
expansion and compression spaces as well as the shaft power output 
corresponding to different operating speeds. Cheng and Yang [12] 
also developed a numerical model to predict the transient varia¬ 
tions of temperatures, pressures and working fluid masses in the 
individual working spaces, so as to obtain the thermodynamic 
behavior of a Stirling engine. Rogdakis et al. [13] analyzed the ther¬ 
modynamic performance of the Solo Stirling Engine V161 unit using 
a computer code based on an adiabatic model. Campos et al. [14] 
developed a dimensionless mathematical model, which combines 
fundamental and empirical correlations, and principles of classical 
thermodynamics, mass and heat transfer accounting for variable 
heat transfer coefficients, to simulate the thermodynamic behavior 
of a Stirling engine. These analyses generally fall into the second 
order analysis according to Martini’s [10] classification, in which the 
engine is divided into several (usually five) chambers. Zero- or one¬ 
dimensional equations of mass continuity, momentum and energy 
conservation are solved. Correlations of various flow friction and 
heat transfer losses are treated as mutually independent and are 
used to modify the total power output. However, in order to gain a 
deeper insight into the complex fluid flow and heat transfer pro¬ 
cesses that occur in the internal gas circuit, and to more accurately 
predict the engine’s performance, a three-dimensional (3D) dy¬ 
namic mesh computational fluid dynamics (CFD) method is rec¬ 
ommended to aid the engine design [15] since it can accommodate 
complex geometries, complicated boundary conditions and variable 
physical properties, depends less on empirical correlations that are 
usually obtained under certain limited experimental conditions. 
With the rapid development of computing power and CFD tools, the 
utilization of CFD method in the research and development of Stir¬ 
ling engines [16], thermoacoustic engines [17], pulse tube re¬ 
frigerators [18] and Stirling regenerator [19] is in an increasing trend 
in recent years. This method is especially useful in this work since no 
exact experimental correlations are readily available for the recip¬ 
rocating flow and heat transfer problem in a hexagonal channel with 
relatively thick wall. 

In this work, a high performance model Stirling engine installed 
with a compact porous-sheets regenerator that has won the first 
prize in the 14th Stirling Techno-rally in Japan is identified for 
theoretical analysis, in order to understand the inherent physical 
mechanism, and to provide guidance for the future scale-up design. 
A commercial CFD code FLUENT is utilized to simulate the recip¬ 
rocating flow and heat transfer in the porous-sheets heat exchanger 


by a dynamic mesh method, and validated with analytical solutions 
and experimental results. Then the flow and heat transfer in the 
entire Stirling engine are numerically simulated using the 3D dy¬ 
namic mesh CFD method. The output power and the thermal effi¬ 
ciency are also obtained. Finally an optimization method is evolved 
by further taking the total entropy generation into account, and 
comparison of comprehensive performance is made between the 
optimized porous-sheets regenerator and the optimized wire mesh 
regenerator. 

2. Description of the model Stirling engine installed with a 
porous-sheets regenerator 

A schematic of the model Stirling engine is shown in Fig. 1(a), 
the detailed drawings of which can be found at Mr. Fukui's website 
[20], The engine is composed of a cold cylinder, a heat exchanger 
and a hot cylinder arranged in a type configuration. The heat 
exchanger includes a heating section, a regenerating section and a 
cooling section, all integrated into the same unit. Two stainless steel 
end caps form two chambers, connecting the heat exchanger with 
the cold cylinder and the hot cylinder. The heating section of the 
heat exchanger together with the hot end cap is heated by a gas 
burner, and the cooling section together with the cold end cap is 
cooled by a heat pipe connected to a heat sink. A Ross drive 
mechanism ensures a phase shift of 90° between the two pistons. 
The bore and stroke of both cylinders are 32.5 mm and 21.4 mm. 
The heat exchanger is constituted by 365 pieces of 0.2-mm-thick, 
circular shaped, porous, stainless steel sheets, each porous sheet 
having 685 hexagonal shaped holes arranged as shown in Fig. 1(b), 
each hole having a side length (a) of 0.4 mm. The porous sheets are 
laminated and inserted into a cylindrical container with an inner 
diameter (4>) of 28.8 mm, the holes of all sheets being aligned with 
each other so as to form 685 flow channels for the working gas, 
each channel having a hexagonal cross section and a length (I) of 
73 mm. The engine is charged with helium at atmospheric pressure. 
The measured rotational speed ( n ) is = 2600 rpm under load when 
the engine is driving a model car. 

3. The reciprocating laminar flow and heat transfer in a 
single channel 

In order to validate the dynamic mesh method used in the CFD 
simulation, a single channel model is created as shown in Fig. 2. The 
channel outer wall is taken along the center planes of the solid 
skeleton that divide the adjacent channels. The calculation domain 
includes a single fluid channel with exactly the same shape and size 
as the actual fluid channel, a solid wall with half the thickness of the 
real sold material and two cylinders with diameters 3 times that of 
the hydraulic diameter of the fluid channel. 

Some general assumptions are made as follows according to the 
practical condition, 

(1) The working fluid is an ideal gas; 

(2) No leakage of working fluid exists; 

(3) Adiabatic boundary condition is specified on the outer wall of 
the engine except at the heating and cooling portions of the 
heat exchanger; 

(4) Sinusoidal movement of the two pistons are specified; 

(5) Cyclic steady state with constant frequency is assumed. 


3.1. The reciprocating laminar flow 

The basic equations for transient fluid flow and heat transfer can 
be found in the user manual of the FLUENT 14.0 software [21], 
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Fig. 1. Schematics of the model Stirling engine and the porous-sheets heat exchanger [ 20], (a) The model Stirling engine, (b) One of the porous sheets. 1 - cold cylinder; 2 — cooling 
end chamber; 3 - porous-sheets heat exchanger; 4 - heating end chamber; 5 - hot cylinder; 6 - Ross drive mechanism. 


which is omitted in this text. In order to validate the dynamic mesh 
CFD model by comparing with available analytical solutions for 
reciprocating laminar pipe flow, the flow field of working gas 
without compression and expansion, without heating and cooling, 
is calculated first. Therefore, the energy equation is deactivated in 
the flow field simulation, and constant thermophysical properties 
at ambient temperature and atmospheric pressure are set for the 
working gas. 

Zhao [22] has presented a laminar—turbulent transition criteria, 
/Seri, for the reciprocating pipe flow as 



hydraulic diameter, m is the oscillatory frequency, Re w = tod^/v is 
the kinetic Reynolds number, where v is the kinematic viscosity of 
fluid. Since the parameter f) of the reciprocating channel flow in 
question is only 93.0, far below the critical value /J cr j of 761 pre¬ 
sented by Zhao [22], the flow is treated as laminar in the 
calculation. 

A “layering” scheme is used for the dynamic mesh method, in 
which the cell quantity is varying and the cell shape is deforming 
with time in both cylinders. A “rigid body” motion is specified to 
each of the two piston walls with the following velocity profile by 
using a compiled User Defined Function (UDF). 

U = U max sin OJt (2) 


where A 0 = x max /d h is the dimensionless oscillation amplitude of 
fluid, where x max is the amplitude of fluid displacement, dh is the 



(a) 


Cooling Regenerating Heating 
Cold section section section 


Hot 



(b) 


Fig. 2. Schematics of the single channel model for validation of the dynamic mesh CFD 
method, (a) Mesh in the channel cross section, (b) A single channel with two pistons. 


where U and U ma x are the transient speed and the maximum piston 
speed. U max is determined by multiplying the real piston speed 
with the cross-sectional area ratio of the piston to the flow channel, 
to ensure that the flow velocity in the channel equals to the real 
flow velocity. All the other boundaries are set as hydrodynamically 
no-slip and thermally adiabatic walls. An initial gauge pressure of 
0 Pa is applied to the whole flow field. 

The PISO scheme is used for the pressure — velocity coupling 
because of its high efficiency in transient flow simulation. A con¬ 
stant time step corresponding to a crank angle of 1° is adopted 
during the simulation. A residual of 1 x 1CT 4 for the continuity 
equation and the momentum equations in 3 directions is set as the 
convergence criteria. The cross-sectional average pressure at the 
channel inlet and outlet and the area-weighted average shear stress 
at the channel inner wall are monitored during the simulation. 

A grid independence test shows that when the initial mesh is 
created with 1.7 million cells, the resulting pressure — time curve 
varies less than 0.1% if further increasing the initial cell quantity. 

After about 10 cycles of iteration, the pressure amplitude varies 
less than 0.1% relative to that of the previous cycle, so the cyclic 
steady state can be considered as having been reached. Then the 
calculation results of flow friction loss and total pressure drop at the 
cyclic steady state are compared with the following analytical solu¬ 
tion for fully developed laminar reciprocating flow in a circular pipe 
derived by Uchida [23] and experimentally validated by Zhao [24], 


= ^P si n(0 + <M (3) 

where <j> is the crank angle, and fa is the phase shift (in degrees) 
between the cross-sectional mean velocity and the wall shear stress, 
Cf,oo is the frictional coefficient for fully developed flow defined by 
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(a) (b) 



(c) 

Fig. 3. Validation of the dynamic mesh model, (a) Re w = 1.07; (b)Re w = 2.14; (c) Re w = 71.9. Num. - numerical; Ana. - analytical; Ap - total pressure drop; 4i w L/d h — frictional loss 
component; pL(du m /dt) — temporal acceleration component. 


m _ Mt) _ -KvV), 

^f, OO v ~/ 1 7 1 o 

lP<ax 2P u max 


(4) 


Ci 


berabei' a - beiaber'a 
ber 2 a + bei 2 a 


(7a) 


where t is time, r w is time varying area weighted average shear 
stress at the channel inner wall, p is the fluid density, u max is the 
maximum cross-sectional mean velocity of the cyclic flow in the 
channel, p is the dynamic viscosity of the fluid, V is the velocity 
vector, F u is calculated by 


yc 2 +c| 

16 v /(a-2C 1 ) 2 + 4C2 
01 is calculated by 


(5) 


^ =tan_ 1 (^ 1 )- tan_ 1 (t) (6) 

where a is the Womersley number defined by a = 0Re w /2, and C\, 
C 2 are constants given by 


_ beraber'a + beiabei'a , 

C 2 = -n- T~2 - 

ber a + bei a 

where ber'a = d(bera)/da and bei'a = d(beio:)/da. 

The comparison of wall shear stress and total pressure drop with 
analytical solutions under different Re u are shown in Fig. 3, in 
which the analytical total pressure drop is calculated by 

Ap = (8) 

Good agreement is found between the analytical solution for 
fully developed flow and the dynamic CFD results with a maximum 
deviation of 2.37%, which can be attributed to the entrance effects 
when the fluid enters and exits the channel, so the analytical so¬ 
lution is to be integrated into the CFD simulation of the overall 
engine in the later sections to reduce the computational cost. The 
parameters and versus Re w in the analytical expression are 
plotted in Fig. 4(a), and the cyclically averaged frictional coefficient 
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Fig. 4. Analytical solution, (a) Parameters F„ and (in degree) in the analytical expression, (b) Cf >0Oi ave- 


is plotted in Fig. 4(b) for the ease of use in the investigation of the 
model Stirling engine, where Cf co.ave is given by 


Q 


, oo ,ave 


64 

nA 0 


( 9 ) 


A residual of 1 x 10~ 4 for the continuity equation and the mo¬ 
mentum equations in 3 directions, and a residual of 1 x 1CT 7 for the 
energy equation are set as the convergence criteria. The other 
solver settings are the same with the previous setting for fluid flow 
simulation. The local convective heat transfer coefficients are ob¬ 
tained by data post-processing for a specific time step using the 
following definition. 


3.2. The conjugate heat transfer with limited temperature 
difference 

Still referring to Fig. 2, in order to evaluate the conjugate heat 
transfer characteristics of the heat exchanger, the heat conduction 
in the solid substrate and the convective heat transfer in the 
reciprocating flow are simultaneously solved. Two types of thermal 
boundary conditions are applied to the outer wall of the heating 
section: 

(1) Constant heat flux of 5000 W/m 2 , corresponding to such 
cases as heating by concentrated solar irradiation; 

(2) Constant temperature of 923 K, corresponding to such cases 
as heating by heat pipe. 


h x 


<?wiM 

T w i(x) - T b (x) 


(11a) 


where q w fx) is the local average heat flux through the inner 
channel wall calculated by 

<?wi(*) =pj <Jwi(x,y,z)dZ (lib) 

r 

where 7’ is the hexagonal boundary curve of the channel cross 
section, and P is the corresponding curve perimeter; T wi (x) is the 
local average temperature of the inner channel wall calculated by 


In both cases, a constant temperature of 350 K is set at the outer 
wall of the cooling section, and an adiabatic boundary condition is 
applied to the outer wall of the regenerating section. All the inner 
walls of the 3 channel sections are set as “Coupled walls”, i.e., hy- 
drodynamically non-slip wall satisfying the temperature continuity 
and heat flux continuity conditions at both sides of the wall, which 
are described as 


7wi(x) = IJ T wi {x,y,z)dl (11c) 

r 

and T b (x) is the local mixed mean temperature defined by 
j pc p T(x,y,z)u(x,y,z)dA J T(x,y,z)u(x,y,z)dA 


(Ts)w ~ 1 

w 

'w 

(10a) 

J pc p u(x.y,z)dA 

J u(x,y,z)dA 

(<Js)w = 1 

M 

>w 

(10b) 

Ac 

A c 

(lid) 


The energy equation is activated, and the idea gas law is 
adopted for the helium density. The thermal conductivity and the 
dynamic viscosity of the fluid, the thermal conductivity and the 
specific heat capacity of the stainless steel substrate are all set as 
piecewise linear functions of temperature. The dynamic mesh and 
the moving piston wall settings are the same as the previous 
settings for the flow field simulation. All the other boundaries are 
set as hydrodynamically non-slip and thermally adiabatic walls. 
An initial condition of zero gauge pressure is applied to the whole 
flow field. 


where p and c p are considered as functions of longitudinal location 
only, and are treated as constant within a local cross section since 
the temperature variation within a cross section is negligible 
compared to the longitudinal variation. 

The local Nusselt number is defined by 


Nu x 


Mh 

k { (x) 


( 12 ) 


where kf(x) is the local thermal conductivity of the working gas. 
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(a) 


(b) 




(c) (d) 

Fig. 5. Heat transfer characteristics of the porous-sheets regenerator under different boundary conditions at the heater wall, (a) h x (0), constant temperature; (b) Nu x (0), constant 
temperature; (c) h x (0), constant heat flux; (d) Nu x (0), constant heat flux. 


The spatial and temporal distributions of h x and Nu x in the 
regenerator section under the two types of thermal boundary 
condition are calculated and illustrated in Fig. 5. It is interesting to 
observe that in both cases, although the local heat transfer coeffi¬ 
cient varies largely and almost linearly in the longitudinal direction 
due to the large difference of gas thermal conductivity at different 
temperature, the local Nusselt number varies not much (less than 
3.7%) and shows a valley-shaped pattern in the longitudinal di¬ 
rection, which might be caused by the influence of entrance and 
exit effects at both ends. It is also observed that both the heat 
transfer coefficients and the Nusselt number varies little (less than 
1.4%) with crank angles except at some crank angles when the fluid 
velocity approaches zero and the data uncertainty might be 
increased, for examples, at 5° and 175°. There is insignificant dif¬ 
ference between the two types of boundary conditions when the 
convective heat transfer coefficient and the Nusselt number are 
concerned, which may indicate that the influence of different 
thermal boundary conditions at the heating section is not easy to be 
transferred to the regenerator section. It is also interesting to find 
that the calculated Nusselt number in both cases is somewhat 
higher than the analytical solutions for fully developed unidirec¬ 
tional laminar flow in a long pipe with hexagonal shaped cross 
section, which has a value of 4.00 for the uniform wall heat flux 
boundary condition, and 3.34 for the uniform wall temperature 
boundary condition [25], This might be explained by the difference 
in the flow and thermal conditions. The reciprocating flow might 


disturb the boundary layer and reduce its thickness, which is 
favorable for enhancing the convective heat transfer. On the other 
hand, the adiabatic boundary condition at the outer channel wall is 
different with either the constant wall heat flux boundary condi¬ 
tion or the constant wall temperature boundary condition, and the 
limited heat capacity of the solid substrate material might be a 
major limit for the heat transfer. 

Based on the calculated results of heat transfer characteristics as 
shown above, a space-cycle averaged Nusselt number defined as 
follows might be appropriate for the engine design in most cases. 

2tu L 

Nu ave = J J Nu x (0)dxd0 (13) 

o o 

In this paper, the calculated Nu ave values are 5.33 for the case of 
constant heater wall heat flux, and 5.31 for the case of constant 
heater wall temperature, respectively. 

4. Dynamic mesh CFD simulation of the entire engine 

4.1. Validation of the porous media sub-model 

Complete meshing of the solid substrate and fluid in the heat 
exchanger portion is not suitable for engineering application or 
even impossible due to the extremely high computational cost. In 
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Fig. 6. Validation of the porous media model, (a) The test calculation model used for validating the porous media model; (b) At 2600 rpm, 1 atm, 300 K 
= 1.07, 1 = 1.2427 x 10 s Jj); (c) At 3500 rpm, 50 atm, 300 K (Re„ = 71.9, 1 = 1.689 x lO 8 ^) 


order to significantly increase the computational efficiency for the 
full engine model, the porous media sub-model is used to integrate 
the analytical solutions of reciprocating laminar flow into the heat 
exchanger portion of the Stirling engine. The determination of the 
parameters in the porous media sub-model is as follows. 

Let 


U m — U m ax sin (f> 


(14) 


where u m is the cyclically averaged cross-sectional mean fluid ve¬ 
locity in the exchanger channel. 

The analytical solution for the frictional loss per unit length in 
the channels is given by 


-4 jr = ~P^P U max ^ + ft) (15) 

u h A max 

The momentum sink term in the porous media model [21] is 
described by 


Sj — — ^ (eu m ) + C 3 — p(t-u m ) 2 

Compare the above two formulae, let 


(16) 


2 

"r“~“P u maxt; / sin(0 + </> 1 )d(</> + ‘/’i) = ^ (eu m ) (17a) 

A max a; j a 


C 3 = 0 (17b) 

where the regenerator porosity e is calculated by e = nA ch /A re g, 
where n is the total number of channels, A C h is the cross-sectional 
area of a single channel, and A reg is the frontal area of the regen¬ 
erator. Then the viscous friction coefficient is obtained as 


1 32 wpFfj 


A test calculation model as shown in Fig. 6(a) including a porous 
media portion and two cylinder portions is designed to validate the 
porous media model. The porous media portion has the same 
dimension with the real heat exchanger, and the frictional coeffi¬ 
cient and porosity parameters are determined according to the 
above method. A sinusoidal movement is specified on each of the 
two pistons without phase shift. The other settings are similar to 
that of the previous sections. The simulation results at ambient 
temperature under low and high kinetic Reynolds numbers are 
shown in Fig. 6(b) and (c). At the low Re w of 1.07, the deviation of the 
porous media model result from the analytical solution is less than 
0.6%, while with the increasing of Re w , the difference between the 
porous media model result and the analytical solution increases, 
reaching 2.9% at a Re w of 71.9, which might be due to the simplifi¬ 
cation by neglecting the inertial loss term as expressed in Eq. (17b). 
Since the Re w value in the real model Stirling engine is close to 1.0, 
the neglect of inertial loss term is acceptable in the porous media 
model. 

In order to validate the dynamic CFD method for the full engine 
model, experiment is conducted on the real model Stirling engine to 
measure the transient pressure drop through the porous-sheets heat 
exchanger. The phase shift between the two pistons is adjusted to 
180° to exclude the effect of compression and expansion. As shown in 
Fig. 7(a), two semiconductor pressure transducers (JTEKT PMS-5M-2- 
50K) are inserted into the two end chambers of the porous-sheets 
heat exchanger. The signal is recorded by a data acquisition system 
(Yokogawa WE7000). The phase angle and rotating frequency is 
measured by an angular encoder, which is controlled by a digital 
signal processor. For each frequency value, after waiting for 30—60 s 
to reach a steady state, data are sampled and stored for 20 revolu¬ 
tions. Air is selected as the working fluid in the experiment and the 
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Fig. 7. Validation of the CFD model with experiment, (a) Experimental setup; (b) Comparison of pressure drop between numerical result and experimental data. 


working frequency ranges from 6.3 to 63 Hz. The pressure drop signal 
is fitted to a sinusoidal form as 

Ap = (Ap) max sin (ip + <j> 2 ) (19) 

Then numerical simulation is conducted using the dynamic 
mesh CFD method incorporating the porous media model. The 
geometric parameters, working fluid and rotating frequency are 


assumed all the same as measured values. The pressure amplitude 
(Ap) max at the corresponding position of the pressure transducer is 
extracted and plotted against rotating frequency in Fig. 7(b) in 
comparison with the experimental data. The experimental error 
mainly comes from the inaccuracy of crank angle and pressure 
measurement. The uncertainty of the crank angle is estimated to be 
±0.25° from the accuracy limitation of the origin setting and the 



Fig. 8. Velocity vectors, (a) $ = 0°; (b) 0 = 90°; (c) 0 = 180°; (d) 0 = 270°. 
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sensitivity of the rotary encoder. The inaccuracy of the pressure 
transducer is less than ±0.3% of the pressure amplitude. The 
simulation results agree quite well with the experimental data with 
a maximum relative error of 4.6%. 

4.2. Flow with compression and expansion under ambient 
temperature 

The computational domain for the full engine model includes 
two cylinders, the heat exchanger, the cooling end chamber that 
connects the cold cylinder with the heat exchanger, and the hot end 
chamber that connects the hot cylinder with the heat exchanger. The 
following motion speed functions are specified for the cold side 
piston and the hot side piston by UDF to take into account the 
expansion and compression effects in the real engine. 

tfcp = Umax sin tot (20a) 

u hp = Umax sin(wt + £) (20b) 

A floating operating pressure is used to reduce the trunca¬ 
tion error in such a closed system with compression and 
expansion. Ideal gas law is used for the gas density, but the 
temperature of all cell zones is fixed to a constant value of 
300 K in the cold state simulation in order to study the effects 
of compression and expansion. The settings for other thermo¬ 
physical properties, boundary conditions and the solver set¬ 
tings are similar to that of the dynamic mesh simulations in 
the previous sections. 

Fig. 8 shows the 3D velocity vectors at different crank angles. 
It can be observed that the flow distribution in the heat 
exchanger is quite uniform due to relatively low pressure drop in 
the connecting chambers and cylinders compared with that in 
the heat exchanger. 

Fig. 9 shows the pressure variation at different locations of the 
engine with time when a cyclic steady state is reached. The pres¬ 
sure difference among different locations can be attributed to the 
internal flow resistance of the gas passage. The simulation result of 
pressure magnitude is compared with the ideal gas law, and the 
deviation is within 1.2%. 

4.3. Flow with compression, expansion, heating and cooling 

To better characterize the heat transfer between the working 
gas and the solid substrate with limited temperature difference and 




Fig. 9. Pressure variation at different locations versus crank angle in a cycle with 
compression and expansion at the constant temperature of 300 K. 


limited heat capacity, a non-thermo-equilibrium model is used for 
the porous media portion. This model is based on a “dual cell” 
approach, which involves a second solid cell zone that overlaps, i.e„ 
spatially coincident with, the porous fluid zone. The two zones are 
solved simultaneously and this solid zone only interacts with the 
fluid with regard to heat transfer. The conservation equation solved 
for the fluid and solid zones can be found in the FLUENT user 
manual [21], which is omitted here. A constant Nusselt number of 
5.3 and a specific heat exchange area of 2.6699 x 10 3 m 2 /m 3 are 
adopted for the non-thermo-equilibrium porous media model. 

The flame is concentrated upon the heating section of the heat 
exchanger assembly, and the combustion gas is directed to the hot 
cylinder, so the hot end of heat exchanger assembly, the hot cyl¬ 
inder and the connecting chamber between them are all exposed to 
hot combustion gas environment, then the combustion gas func¬ 
tions as fine thermal insulator to prevent heat loss, therefore in the 
CFD modeling it is assumed that heat enters and leaves the system 
only through the heating and cooling sections of the heat 
exchanger, while all the other boundaries are treated as adiabatic 
walls. The outer surface of the heating portion is assigned a con¬ 
stant heating power, Qj n , of 76.0 W, and the outer surface of the 
cooling portion is assigned a constant temperature of 330 K as the 
thermal boundary conditions. The initial temperature field is 



(a) 


(b) 


Fig. 10. Pressure variation with compression and expansion with heating and cooling, (a) Pressure variation at different locations versus crank angle in a cycle, (b) p-V diagram. 
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Table 1 

Engine specifications and calculated performance data. 


Bore x stroke (mm) 


Gas mean temperature 

(K) 


Expansion cylinder 

32.5 x 21.4 

Heating end 

367 

Compression cylinder 

32.5 x 21.4 

Cooling end 

1020 

Dead volume (cm 3 ) 

25.50 

Rotating speed 
(rpm) 

2600 

Heating end chamber 

1.72 

Pressure ratio 
(Pmax/Pmin) 

2.15 

Heat exchanger 

20.78 

Compression ratio 
(V max /V mi „) 

1.64 

Cooling end chamber 

1.72 

Power output (W) 

33.4 

Clearance and conduit 

1.28 

Specific power 
(W/cm 3 ) 

1.88 

Mean pressure P m (bar) 

1.48 

Indicated thermal 
efficiency (%) 

43.9 


obtained by solving a steady state energy equation before the 
transient calculation begins, and applied to the overall calculation 
domain as an initial thermal condition. The initial absolute pressure 
is set to a value of 1.3 x 10 5 Pa. The settings for the thermophysical 
properties, the dynamic mesh, the other boundary conditions and 
the solver are similar to that of Section 4.2. 

Fig. 10(a) shows the pressure variation at different locations of the 
working space against time when a cyclic steady state is reached. The 
pressure difference among different locations can also be attributed 
to the internal flow resistance of the gas passage. The area-weighted 
pressure of working gas adjacent to the two pistons is plotted against 
the variation of total volume, and is shown in Fig. 10(b). The main 
specifications and calculated performance of the model Stirling en¬ 
gine are shown in Table 1. An indicated power of 33.4 W is obtained 
by integrating the p—V diagram, corresponding to a specific power of 
1.88 W/cm 3 based on the swept volume of the expansion cylinder, 
and an indicated thermal efficiency (?) = W/Qj n ) of 43.9%, which is 
considered high performance among model Stirling engines under 
low operating pressure. 


diameter <F, calculate the thermal properties p s , pf, c Pi f, c s , k s , 
kf, p, the mean mass flow rate m, the mean fluid temperature 
Tf m and the blow time tbiow for half cycle. 

(2) Assume a mesh wire diameter, d w , then calculate Nusselt 
number, Nu, Biot number, Bi and Fourier number, Fo. 

Nu d = 0.42Re d 56 , for mesh wire [26] (21a) 

JVu dh = 5.3, for hexagonal porous sheets (21b) 

where Nu d = hd w /k f , Re d = p f u m d w /p, where u m is the flow ve¬ 
locity within the matrix defined by u m = u 0 /fl, where u 0 is the 
frontal flow velocity, and /? = A f /A reg is ratio of the free flow area Af 
to the frontal area A re g- The Biot number is defined as 
Bi = hd w /2k w for mesh wire, and Bi = h5 w /k w for porous sheets. 
The Fourier number is Fo = 4a w t blow /d 3 , for mesh wire, and Fo = 
Gwfbiow/^w f° r porous sheets. 

(3) Select an appropriate d m , to ensure 0.95 < 

0(Bi. Fo)| x=0 or r=0 < 1 according to the following analytical 
solutions [28], where 0 = (T - T 0 )/(T [m - T 0 ) is the 
dimensionless centerline temperature of the solid material 
(mesh wire or porous sheet), T is the transient local solid 
temperature, Tf m is the mean fluid temperature and To is the 
initial solid temperature. 


@ lx=0, or r=0 = 1 - I] C n exp(~MnF°) 
n = 1 

For porous sheets/parallel plates, 


C n 


2 sin p n 

p n + cos p n sin p n 


(22a) 


(22b) 


5. Optimization of porous-sheets regenerator and 
comparison with optimized wire mesh regenerator 

In this section, the authors evolve the regenerator optimization 
method originally presented by Hamaguchi etc. [26], which had 
been validated with real engine tests by Miyabe et al. [27], to 
further incorporate the total entropy generation rate into consid¬ 
eration. Given the same operating condition, the geometrical pa¬ 
rameters of both porous sheets regenerator and conventional wire 
mesh regenerator can be optimized in terms of heat transfer 
characteristics, flow resistance and entropy generation rate, so the 
comprehensive performance of the two types of regenerators can 
be compared on an equal basis. The optimization procedure is 
briefly described as follows. 


where p n is positive roots of the following transcendental equation. 


p n tan p n = Bi, n = 1,2, 
For cylindrical mesh wires, 

Cn = 2Ji (p n ) 

Mn[io(f i n)+Jl(Mn)] 


(22c) 


(22d) 


where p n is positive roots of the following transcendental equation. 


FnJ l(pn) ^ — Ij 2, ••• 


(22e) 


(1) Given the same combination of operating conditions, 
including the mean pressure P m , the hot and cold end tem¬ 
perature Th and T c , the rotational velocity n and the matrix 


(4) Select number of layers for wire mesh or porous sheets, N to 
satisfy e(N TU , Cr) > 0.95 according to the following analytical 
solution [26], 


N TU 




Wlo(^) + j (2 + 'g!- 2 ) t ->"»->/ 0 (2vtw)d7 


e 


(23) 
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Table 2 

The geometric parameters of optimized wire mesh [27], 


Mesh # 

d w [mm] 

p [mm] 

/ [mm] 

13 

<P 

a [mm 2 /mm 3 ] 

300 

0.040 

0.084 

0.044 

0.274 

0.586 

32.52 

200 

0.050 

0.127 

0.077 

0.368 

0.668 

21.72 

150 

0.071 

0.169 

0.098 

0.336 

0.642 

16.26 


Note: p — pitch; / — aperture side length; a — specific surface area. 


Table 3 

The geometric parameters of optimized porous sheets. Sheet thickness 5 t = 0.1 mm. 


Porous sheet 

2<5 W 

a [mm] 

■©- 

II 

'Q. 

a = —— -j [mm 2 /mm 3 ] 

code # 

[mm] 


a 2 

73 ( 0 +^ 




("+A 5 ") 


A 

0.022 

0.060 

0.681 

26.22 

B 

0.048 

0.075 

0.533 

16.42 

C 

0.058 

0.100 

0.561 

12.96 

D 

0.024 

0.100 

0.771 

17.81 


Note: a — side length of the hexagonal cross section; <5 W - wall thickness as denoted 
in Fig. 1(b). 


where e is the thermal effectiveness of the regenerator, 
e = (Tp, - T fc )/(T mh - T fc ), where T ni and T fc are hot and cold end 
fluid temperature, and T m h is the hot end matrix temperature, N TU 
is the number of heat transfer units, N TU = ( hS)/(mc p ), where S is 
the total heat transfer area between solid matrix and fluid, m is 
the mean mass flow rate. Cr = MC m /(mc p f b | 0W ) is the heat 
capacity ratio of the solid matrix to the fluid flowing through the 
matrix. 


(5) Readjust N to minimize the total entropy generation rate, 
S g = S g h + S g f, where S g h = Q\.^_ - 7 jl ~Q.fi- is the portion 
caused by the irreversible heat transfer between solid matrix 
and working fluid, where Q is the total heat transfer between 
them, and AT = |.T f — T s |.; S gf is caused by the flow friction. 
Also check that the flow frictional loss (Ap) f is within 
acceptable range. 


AT = 


For a porous-sheets regenerator, 

Q d h 


SNu db k f 
1 


-*g,f ~ 2 00 ,ave^lRm| 

, 4L 1 2 

(Ap)f — Q, oo,ave2 P u max 


(24a) 

(24b) 

(24c) 


where S is the total heat transfer area between matrix and fluid. 
For a wire mesh regenerator, 


AT =4 d ' 


S Nu d kf 


(25a) 


% = (Ap) f A 0 |Lfo|^ 


(25b) 


(Ap) f = N^pulJ (25c) 

where N is the number of mesh screens, and / is the flow friction 
coefficient given by the experimental correlation [27] 

/ = ^+0.337 (25d) 

where Re ( is defined by Re ( = pu m //p, where / is the side length of 
the mesh aperture. 

Some typical operating conditions are selected for optimization 
according to the above procedure. The resulting optimum geometric 
parameters of wire mesh regenerators and porous-sheets re¬ 
generators are shown in Tables 2 and 3. The major operating param¬ 
eters are compared in Table 4. It is found that under each operating 
condition, the porous-sheets regenerator has significantly lower flow 
frictional loss while keeping the same or higher thermal effectiveness, 
thus leading to significantly less total entropy generation rate (38%— 
51% less) in comparison with the corresponding wire mesh regener¬ 
ator. Since the total entropy generation rate, S g is an indicator of 
available work loss, W ( = T 0 S g , where To is the absolute ambient 
temperature, it can be reasonably used to measure the work loss 
within the regenerator caused by irreversible heat transfer and flow 
friction, therefore it should be minimized in the engine design in order 
to obtain maximum power output and thermal efficiency. 

The step (5) of the above optimization procedure is based on 
minimization of total entropy generation rate, then the resulted 
dead volume may be different for the two kinds of regenerator. If 
the number of matrix layers, N is selected to let the two kinds of 
regenerator have equal dead volumes while assuming all the same 
operating conditions as in Table 4, the resulted thermal effective¬ 
ness, pressure drop and total entropy generation rates are 
compared in Table 5. It can be found that the total entropy gener¬ 
ation rates of the porous-sheets regenerator only amount to 34.8%— 
65.8% of that of the corresponding wire mesh regenerator. In other 
words, the frictional pressure drop through the porous-sheets 


Table 4 

Comparison of major parameters between optimized porous-sheets regenerator and optimized wire mesh regenerator (T c = 400 K, T e = 900 I<). 


Operating condition 

P m = 1.879 MPa 

P m = 1.957 MPa 

P m = 3.132 MPa 

P m = 1.253 MPa 

n - 2600 rpm 


n = 1000 rpm 


n = 600 rpm 


n = 1200 rpm 


0 = 28.8 mm 


0 = 28.8 mm 


0 = 28.8 mm 


0 = 14.0 mm 


Regenerator element 

Wire mesh 

Porous sheet 

Wire mesh 

Porous sheet 

Wire mesh 

Porous sheet 

Wire mesh 

Porous sheet 

M300 

A 

M200 

B 

Ml 50 

C 

M150 

D 

Bi 

0.0117 

0.00709 

0.00672 

0.0124 

0.00841 

0.0112 

0.0167 

0.00387 

Fo 

134.6 

444.9 

223.9 

243.0 

185.1 

277.3 

92.5 

809.9 

Tr 

0.956 

0.957 

0.950 

0.950 

0.955 

0.955 

0.953 

0.956 

Cr 

80.4 

141.0 

206.0 

323.2 

347.3 

540.4 

52.3 

72.8 

Ntu 

200.0 

405.5 

511.2 

826.9 

878.8 

1451.4 

131.0 

214.9 

N 

107 

195 

285 

318 

502 

905 

128 

396 

e 

0.982 

0.990 

0.993 

0.997 

0.998 

0.999 

0.973 

0.982 

(Ap)f(kPa) 

13.86 

6.822 

5.651 

3.498 

5.269 

3.192 

14.03 

8.601 

S g (W/K) 

0.0657 

0.0324 

0.0103 

0.0064 

0.0058 

0.0035 

0.0308 

0.0188 


Note: P m — mean operating pressure; n — rotational speed; <P — matrix diameter. 
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Table 5 

Comparison between porous-sheets regenerator and wire mesh regenerator based on equal dead volume (T c = 400 K, T e = 900 K). 


Operating 

condition 

P m = 1.879 MPa 

P m = 1.957 MPa 

P m = 3.132 MPa 

P m = 1.253 MPa 

n = 2600 rpm 

n = 1000 rpm 

n = 600 rpm 


n = 1200 rpm 

0 = 28.8 mm 

0 = 28.8 mm 

0 = 28.8 mm 

0 = 14.0 mm 

Regenerator 

element 

Wire 

mesh 

Porous 

sheet 

Wire 

mesh 

Porous 

sheet 

Wire 

mesh 

Porous 

sheet 

Wire 

mesh 

Porous 

sheet 

M300 

A 

M200 

B 

M150 

C 

M150 

D 

V'd.r (cm 3 ) 

7.45 

8.61 

10.97 

5.98 

V d .t (cm 3 ) 

12.17 

13.33 

15.69 

10.69 

N 

244 

168 

198 

248 

185 

300 

426 

503 

e 

0.990 

0.988 

0.988 

0.993 

0.987 

0.990 

0.989 

0.985 

(Ap) f (kPa) 

31.59 

5.877 

3.926 

2.728 

1.942 

1.058 

46.69 

10.92 

S g (W/K) 

0.0892 

0.0327 

0.0110 

0.0066 

0.0089 

0.0058 

0.0557 

0.0194 


Note: Vij'i - total dead volume; V DR — regenerator dead volume. 


regenerator only amounts to a small fraction (18.6%—69.5%) of that 
of the corresponding wire mesh regenerator while keeping similar 
thermal effectiveness. 


Grant No. 51161140332. The authors would also express their 
gratitude to Mr. Fukui of the Fujifilm Co. Ltd., the designer of the 
porous-sheets regenerator in this research. 


6. Conclusion 

In this study, the reciprocating flow and heat transfer charac¬ 
teristics in a model Stirling engine with a compact porous-sheets 
heat exchanger are numerically simulated by a dynamic mesh 
CFD method and are validated with analytical solutions and 
experimental results. The following major conclusions can be 
drawn from the analysis. 

(1) The flow frictional loss in the porous-sheets regenerator with 
regular shaped flow channels under the working condition of 
1 atm and 2600 rpm can be as low as 800 Pa. 

(2) Although the local heat transfer coefficient varies largely in 
the longitudinal direction due to the strong dependence of 
gas thermal conductivity on temperature, the local Nusselt 
number varies less than 3.7% and shows a valley-shaped 
pattern in the longitudinal direction. The Nusselt number 
varies less than 1.4% with crank angles. A spatial-temporally 
averaged Nusselt number of 5.33 for the case of constant 
heater wall heat flux, and 5.31 for the case of constant heater 
wall temperature are obtained respectively, which are 
somewhat higher than the analytical values of 4.00 and 3.34 
for fully developed unidirectional laminar flow in a long pipe 
with hexagonal shaped cross section. 

(3) A calculated power of 33.4 W, corresponding to a specific 
power of 1.88 W/cm 3 based on the swept volume of the 
expansion cylinder, and an indicated thermal efficiency of 
43.9% is obtained. The high performance of the model Stirling 
engine can be attributed to the low frictional loss and the 
excellent heat transfer characteristics, as well as the compact 
heat exchanger design that significantly reduces the dead 
volume. 

(4) Optimization of both the porous sheets regenerator and the 
conventional wire mesh regenerator shows that under given 
operating conditions, the porous-sheets regenerator has 
38%—51% lower total entropy generation rate while main¬ 
taining the same or higher thermal effectiveness, thus lead¬ 
ing to less available work loss, contributing to higher power 
output and thermal efficiency. 
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